Soil Data Exploration

Opening and exploring data

Open soil and elevation data:

#opening csv
soils <- read_csv(here("data", "shift_soil_data.csv")) %>% # read in soil csv
  mutate(date = lubridate::mdy(date), # change to date class
         week = lubridate::week(date)) %>% # create a week column
  select(soil_id:transect, meter_location, electro_cond_mS_per_cm, season, latitude, longitude)%>% #select certain column to work with
  mutate(paired_flight = lubridate::ymd(paired_flight)) %>% 
  drop_na() %>%  #drop any NA rows
  st_as_sf(coords = c("longitude", "latitude")) %>% 
  st_set_crs(value = 4326) %>% 
  st_transform(crs = 32611)


# Open spatial Data
soils_sf <- read_sf(here("data", "soil_w_elevation.shp")) %>% # read in sf file
  mutate(date = lubridate::ymd(date), #change date to date class
         week = lubridate::week(date)) %>% #create week
  select(electro_co, date, transect, meter_loca) %>% # select certain columns
  filter(date != "2022-09-12") %>% 
  mutate(electro_co = as.numeric(electro_co)) %>% #change to numeric class
  drop_na() #drop NAs

bbox <- read_sf(here("data", "bbox.shp")) #read in bounding box file

boundary <- read_sf(here("data", "boundary.shp")) %>% 
  st_transform(crs = 32611)

dem <- read_stars(here("data", "dem.tif")) %>% # read in dem 
  st_crop(y = st_bbox(bbox)) %>%  #crop to bounding box
  st_warp(cellsize = 4.8, crs = 32611)

mllw <- dem + 0.042

mllw_all <- c(mllw, mllw, mllw, mllw, mllw, mllw, mllw, mllw, mllw, mllw, mllw, mllw, mllw, along = 3) %>%
  st_set_dimensions(3, values = lubridate::ymd(c("2022_02_24", "2022_02_28", "2022_03_08", "2022_03_16", "2022_03_22", "2022_04_05", "2022_04_12", "2022_04_20", "2022_04_29", "2022_05_03", "2022_05_11", "2022_05_17", "2022_05_29")), names = "date")

mllw_extract <- mllw %>% # call dem
  st_extract(soils$geometry) # extract elevation to point layer

Join soil and elevation data:

soils <- cbind(soils, mllw_extract$dem.tif) %>% #bind extracted values to the soil data
  rename("elevation" = "mllw_extract.dem.tif") %>% 
  st_transform(crs = 32611)

Descriptive statistics:

#make descriptive stat table
soils %>% # to soils
  group_by(transect) %>%# group by transect  
  summarize(min = min(electro_cond_mS_per_cm),#calculate mins
            max = max(electro_cond_mS_per_cm),#calculate max
            mean = mean(electro_cond_mS_per_cm),#calculate mean
            dev = sd(electro_cond_mS_per_cm),#obtain standard dev
            variance = var(electro_cond_mS_per_cm)) %>% # calculate variance
  kableExtra::kable() %>% # create a table format
  kableExtra::kable_classic(lightable_options = "striped")#theme the table

transect	min	max	mean	dev	variance	geometry
COPR_1_EW	0.42	1.34	0.6941176	0.2020538	0.0408257	MULTIPOINT ((235669.2 38115…
COPR_2_EW	6.32	10.53	8.2875000	1.2435749	1.5464786	MULTIPOINT ((235855 3812030…
COPR_2_NS	6.30	13.54	9.7581250	2.1347107	4.5569896	MULTIPOINT ((235855 3812030…
NCOS_1_NS	0.12	7.47	2.7527778	2.0644399	4.2619121	MULTIPOINT ((235777.4 38125…
NCOS_2_EW	0.74	8.62	2.9125714	1.8626205	3.4693550	MULTIPOINT ((235671.2 38123…
NCOS_2_NS	2.88	10.14	4.6538889	1.4925439	2.2276873	MULTIPOINT ((235685.7 38123…

# make descriptive stat dataframe
soil_stats <- soils %>% #call soils
  group_by(transect) %>% #group by transects
  summarize(mean = mean(electro_cond_mS_per_cm),# create dataframe with these stats
            max = max(electro_cond_mS_per_cm),
            min = min(electro_cond_mS_per_cm),
            range = max - min)

Plot salinity v elevation

ggplot(soils, aes(x = elevation, y = electro_cond_mS_per_cm, color = transect)) +#ggplot of soils/dem data
  geom_point()+#geometry of the plot
  scale_color_manual(values = calecopal::cal_palette(name = "superbloom3", n =6, type = "discrete"))+# colors to use
  labs(x = "Elevation (m)",#labels
       y = "Electrical conductivity (mS/cm)")+
  theme(legend.position = "none")+#no legend theme
  guides(color = guide_legend(title = "Transect ID"))+#changing legend title
  ggtitle("Soil Electrical Conductivity Across Elevations")+# title of plot
  theme(plot.title = element_text(color = "#5b4f41", size = 20, hjust = 0.5),# plot theming
        plot.background = element_rect("white"),
        panel.background = element_rect("#faf7f2"),
        panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
        axis.text = element_text(color = "#5b4f41", size = 14),
        axis.title = element_text(color = "#5b4f41", size = 16),
        strip.background = element_rect("#f8f8f8"),
        strip.text = element_text(color = "#5b4f41"),
        axis.line = element_line(color = "#5b4f41"),
        legend.text = element_text(color = "#5b4f41"),
        legend.title = element_text(color = "#5b4f41"))

Salinity vs. elevation, facet wrapped

strip_labels <- c("COPR_1_EW" = "COPR 1 EW",#create better formated labels for facet wrap strips
                  "COPR_2_EW" = "COPR 2 EW",
                  "COPR_2_NS" = "COPR 2 NS",
                  "NCOS_1_NS" = "NCOS 1 NS",
                  "NCOS_2_EW" = "NCOS 2 EW",
                  "NCOS_2_NS" = "NCOS 2 NS")

ggplot(soils, aes(x = elevation, y = electro_cond_mS_per_cm, color = transect)) +#start ggplot
  geom_point()+# point geometry
  scale_color_manual(values = calecopal::cal_palette(name = "superbloom3", n =6, type = "discrete"))+#colors
  facet_wrap(~transect, labeller = as_labeller(strip_labels))+#subplots based on transect
  theme(legend.position = "none")+# legend theme
  labs(x = "Elevation (m)",#labels
       y = "Electrical conductivity (mS/cm)")+
  guides(color = guide_legend(title = "Transect ID"))+#legend title
  ggtitle("Soil Electrical Conductivity Across Elevations")+#plot title
  theme(plot.title = element_text(color = "#5b4f41", size = 16, hjust = 0.5),#plot theme
        plot.background = element_rect("white"),
        panel.background = element_rect("#faf7f2"),
        panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
        axis.text = element_text(color = "#5b4f41"),
        axis.title = element_text(color = "#5b4f41"),
        strip.background = element_rect("#f8f8f8"),
        strip.text = element_text(color = "#5b4f41"),
        axis.line = element_line(color = "#5b4f41"),
        legend.text = element_text(color = "#5b4f41"),
        legend.title = element_text(color = "#5b4f41"))

Open Raster Files:

Opening and Compiling NDVI:

ndvi_1 <- read_stars(here("imagery", "NDVI", "2022_02_24_ndvi.tif")) %>% 
  setNames("ndvi")
ndvi_2 <- read_stars(here("imagery", "NDVI", "2022_02_28_ndvi.tif")) %>% 
  setNames("ndvi")
ndvi_3 <- read_stars(here("imagery", "NDVI", "2022_03_08_ndvi.tif")) %>% 
  setNames("ndvi")
ndvi_4 <- read_stars(here("imagery", "NDVI", "2022_03_16_ndvi.tif")) %>% 
  setNames("ndvi")
ndvi_5 <- read_stars(here("imagery", "NDVI", "2022_03_22_ndvi.tif")) %>% 
  setNames("ndvi")
ndvi_6 <- read_stars(here("imagery", "NDVI", "2022_04_05_ndvi.tif")) %>% 
  setNames("ndvi")
ndvi_7 <- read_stars(here("imagery", "NDVI", "2022_04_12_ndvi.tif")) %>% 
  setNames("ndvi")
ndvi_8 <- read_stars(here("imagery", "NDVI", "2022_04_20_ndvi.tif")) %>% 
  setNames("ndvi")
ndvi_9 <- read_stars(here("imagery", "NDVI", "2022_04_29_ndvi.tif")) %>% 
  setNames("ndvi")
ndvi_10 <- read_stars(here("imagery", "NDVI", "2022_05_03_ndvi.tif")) %>% 
  setNames("ndvi")
ndvi_11 <- read_stars(here("imagery", "NDVI", "2022_05_11_ndvi.tif")) %>% 
  setNames("ndvi")
#ndvi_12 <- read_stars(here("imagery", "NDVI", "2022_05_12_ndvi.tif")) %>% 
 # setNames("ndvi")
ndvi_13 <- read_stars(here("imagery", "NDVI", "2022_05_17_ndvi.tif")) %>% 
  setNames("ndvi")
ndvi_14 <- read_stars(here("imagery", "NDVI", "2022_05_29_ndvi.tif")) %>% 
  setNames("ndvi")


ndvi_all <- c(ndvi_1, ndvi_2, ndvi_3, ndvi_4, ndvi_5, ndvi_6, ndvi_7, ndvi_8, ndvi_9, ndvi_10, ndvi_11, ndvi_13, ndvi_14, along = 3) %>% 
  st_set_dimensions(3, values = lubridate::ymd(c("2022_02_24", "2022_02_28", "2022_03_08", "2022_03_16", "2022_03_22", "2022_04_05", "2022_04_12", "2022_04_20", "2022_04_29", "2022_05_03", "2022_05_11", "2022_05_17", "2022_05_29")), names = "date")

Opening and Compiling VSSI:

vssi_1 <- read_stars(here("imagery", "VSSI", "02_24_vssi.tif")) %>% 
  setNames("vssi")
vssi_2 <- read_stars(here("imagery", "VSSI", "02_28_vssi.tif")) %>% 
  setNames("vssi")
vssi_3 <- read_stars(here("imagery", "VSSI", "03_08_vssi.tif")) %>% 
  setNames("vssi")
vssi_4 <- read_stars(here("imagery", "VSSI", "03_16_vssi.tif")) %>% 
  setNames("vssi")
vssi_5 <- read_stars(here("imagery", "VSSI", "03_22_vssi.tif")) %>% 
  setNames("vssi")
vssi_6 <- read_stars(here("imagery", "VSSI", "04_05_vssi.tif")) %>% 
  setNames("vssi")
vssi_7 <- read_stars(here("imagery", "VSSI", "04_12_vssi.tif")) %>% 
  setNames("vssi")
vssi_8 <- read_stars(here("imagery", "VSSI", "04_20_vssi.tif")) %>% 
  setNames("vssi")
vssi_9 <- read_stars(here("imagery", "VSSI", "04_29_vssi.tif")) %>% 
  setNames("vssi")
vssi_10 <- read_stars(here("imagery", "VSSI", "05_03_vssi.tif")) %>% 
  setNames("vssi")
vssi_11 <- read_stars(here("imagery", "VSSI", "05_11_vssi.tif")) %>% 
  setNames("vssi")
#vssi_12 <- read_stars(here("imagery", "VSSI", "05_12_vssi.tif")) %>% 
 # setNames("vssi")
vssi_13 <- read_stars(here("imagery", "VSSI", "05_17_vssi.tif")) %>% 
  setNames("vssi")
vssi_14 <- read_stars(here("imagery", "VSSI", "05_29_vssi")) %>% 
  setNames("vssi")

vssi_all <- c(vssi_1, vssi_2, vssi_3, vssi_4, vssi_5, vssi_6, vssi_7, vssi_8, vssi_9, vssi_10, vssi_11, vssi_13, vssi_14, along = 3) %>%
  st_set_dimensions(3, values = lubridate::ymd(c("2022_02_24", "2022_02_28", "2022_03_08", "2022_03_16", "2022_03_22", "2022_04_05", "2022_04_12", "2022_04_20", "2022_04_29", "2022_05_03", "2022_05_11", "2022_05_17", "2022_05_29")), names = "date")

Opening and Compiling mARI:

mari_1 <- read_stars(here("imagery", "mARI", "02_24_mari_avg.tif")) %>% 
  setNames("mari")
mari_2 <- read_stars(here("imagery", "mARI", "02_28_mari_avg.tif")) %>% 
  setNames("mari")
mari_3 <- read_stars(here("imagery", "mARI", "03_08_mari_avg.tif")) %>% 
  setNames("mari")
mari_4 <- read_stars(here("imagery", "mARI", "03_16_mari_avg.tif")) %>% 
  setNames("mari")
mari_5 <- read_stars(here("imagery", "mARI", "03_22_mari_avg.tif")) %>% 
  setNames("mari")
mari_6 <- read_stars(here("imagery", "mARI", "04_05_mari_avg.tif")) %>% 
  setNames("mari")
mari_7 <- read_stars(here("imagery", "mARI", "04_12_mari_avg.tif")) %>% 
  setNames("mari")
mari_8 <- read_stars(here("imagery", "mARI", "04_20_mari_avg.tif")) %>% 
  setNames("mari")
mari_9 <- read_stars(here("imagery", "mARI", "04_29_mari_avg.tif")) %>% 
  setNames("mari")
mari_10 <- read_stars(here("imagery", "mARI", "05_03_mari_avg.tif")) %>% 
  setNames("mari")
mari_11 <- read_stars(here("imagery", "mARI", "05_11_mari_avg.tif")) %>% 
  setNames("mari")
#mari_12 <- read_stars(here("imagery", "mARI", "05_12_mari_avg.tif")) %>% 
 # setNames("mari")
mari_13 <- read_stars(here("imagery", "mARI", "05_17_mari_avg.tif")) %>% 
  setNames("mari")
mari_14 <- read_stars(here("imagery", "mARI", "05_29_mari_avg.tif")) %>% 
  setNames("mari")

mari_all <- c(mari_1, mari_2, mari_3, mari_4, mari_5, mari_6, mari_7, mari_8, mari_9, mari_10, mari_11, mari_13, mari_14, along = 3) %>%
  st_set_dimensions(3, values = lubridate::ymd(c("2022_02_24", "2022_02_28", "2022_03_08", "2022_03_16", "2022_03_22", "2022_04_05", "2022_04_12", "2022_04_20", "2022_04_29", "2022_05_03", "2022_05_11", "2022_05_17", "2022_05_29")), names = "date")

Opening and Compiling CRSI:

crsi_1 <- read_stars(here("imagery", "CRSI", "02_24_crsi.tif")) %>% 
  setNames("crsi")
crsi_2 <- read_stars(here("imagery", "CRSI", "02_28_crsi.tif")) %>% 
  setNames("crsi")
crsi_3 <- read_stars(here("imagery", "CRSI", "03_08_crsi.tif")) %>% 
  setNames("crsi")
crsi_4 <- read_stars(here("imagery", "CRSI", "03_16_crsi.tif")) %>% 
  setNames("crsi")
crsi_5 <- read_stars(here("imagery", "CRSI", "03_22_crsi.tif")) %>% 
  setNames("crsi")
crsi_6 <- read_stars(here("imagery", "CRSI", "04_05_crsi.tif")) %>% 
  setNames("crsi")
crsi_7 <- read_stars(here("imagery", "CRSI", "04_12_crsi.tif")) %>% 
  setNames("crsi")
crsi_8 <- read_stars(here("imagery", "CRSI", "04_20_crsi.tif")) %>% 
  setNames("crsi")
crsi_9 <- read_stars(here("imagery", "CRSI", "04_29_crsi.tif")) %>% 
  setNames("crsi")
crsi_10 <- read_stars(here("imagery", "CRSI", "05_03_crsi.tif")) %>% 
  setNames("crsi")
crsi_11 <- read_stars(here("imagery", "CRSI", "05_11_crsi.tif")) %>% 
  setNames("crsi")
#crsi_12 <- read_stars(here("imagery", "CRSI", "05_12_crsi.tif")) %>% 
 # setNames("crsi")
crsi_13 <- read_stars(here("imagery", "CRSI", "05_17_crsi.tif")) %>% 
  setNames("crsi")
crsi_14 <- read_stars(here("imagery", "CRSI", "05_29_crsi.tif")) %>% 
  setNames("crsi")

crsi_all <- c(crsi_1, crsi_2, crsi_3, crsi_4, crsi_5, crsi_6, crsi_7, crsi_8, crsi_9, crsi_10, crsi_11, crsi_13, crsi_14, along = 3) %>%
  st_set_dimensions(3, values = lubridate::ymd(c("2022_02_24", "2022_02_28", "2022_03_08", "2022_03_16", "2022_03_22", "2022_04_05", "2022_04_12", "2022_04_20", "2022_04_29", "2022_05_03", "2022_05_11", "2022_05_17", "2022_05_29")), names = "date")

Compiling NIR Bands

nir_1 <- read_stars(here("imagery", "subset_images", "rfl_2022_02_24.tif")) %>% 
  st_set_dimensions(3, values = seq(1:425)) %>% 
  filter(band %in% c(76:126)) %>% 
  split(3) %>% 
  setNames(c(76:126))

nir_2 <- read_stars(here("imagery", "subset_images", "rfl_2022_02_28.tif")) %>%
  st_set_dimensions(3, values = seq(1:425)) %>%
  filter(band %in% c(76:126)) %>%
  split(3) %>%
  setNames(c(76:126))

nir_3 <- read_stars(here("imagery", "subset_images", "rfl_2022_03_08.tif")) %>%
  st_set_dimensions(3, values = seq(1:425)) %>%
  filter(band %in% c(76:126)) %>%
  split(3) %>%
  setNames(c(76:126))

nir_4 <- read_stars(here("imagery", "subset_images", "rfl_2022_03_16.tif")) %>%
  st_set_dimensions(3, values = seq(1:425)) %>%
  filter(band %in% c(76:126)) %>%
  split(3) %>%
  setNames(c(76:126))

nir_5 <- read_stars(here("imagery", "subset_images", "rfl_2022_03_22.tif")) %>%
  st_set_dimensions(3, values = seq(1:425)) %>%
  filter(band %in% c(76:126)) %>%
  split(3) %>%
  setNames(c(76:126))

nir_6 <- read_stars(here("imagery", "subset_images", "rfl_2022_04_05.tif")) %>%
  st_set_dimensions(3, values = seq(1:425)) %>%
  filter(band %in% c(76:126)) %>%
  split(3) %>%
  setNames(c(76:126))

nir_7 <- read_stars(here("imagery", "subset_images", "rfl_2022_04_12.tif")) %>%
  st_set_dimensions(3, values = seq(1:425)) %>%
  filter(band %in% c(76:126)) %>%
  split(3) %>%
  setNames(c(76:126))

nir_8 <- read_stars(here("imagery", "subset_images", "rfl_2022_04_20.tif")) %>%
  st_set_dimensions(3, values = seq(1:425)) %>%
  filter(band %in% c(76:126)) %>%
  split(3) %>%
  setNames(c(76:126))

nir_9 <- read_stars(here("imagery", "subset_images", "rfl_2022_04_29.tif")) %>%
  st_set_dimensions(3, values = seq(1:425)) %>%
  filter(band %in% c(76:126)) %>%
  split(3) %>%
  setNames(c(76:126))

nir_10 <- read_stars(here("imagery", "subset_images", "rfl_2022_05_03.tif")) %>%
  st_set_dimensions(3, values = seq(1:425)) %>%
  filter(band %in% c(76:126)) %>%
  split(3) %>%
  setNames(c(76:126))

nir_11 <- read_stars(here("imagery", "subset_images", "rfl_2022_05_11.tif")) %>%
  st_set_dimensions(3, values = seq(1:425)) %>%
  filter(band %in% c(76:126)) %>%
  split(3) %>%
  setNames(c(76:126))

#nir_12 <- read_stars(here("imagery", "subset_images", "rfl_2022_05_12.tif")) %>%
#  st_set_dimensions(3, values = seq(1:425)) %>%
#  filter(band %in% c(76:126)) %>%
#  split(3) %>%
#  setNames(c(76:126))

nir_13 <- read_stars(here("imagery", "subset_images", "rfl_2022_05_17.tif")) %>%
  st_set_dimensions(3, values = seq(1:425)) %>%
  filter(band %in% c(76:126)) %>%
  split(3) %>%
  setNames(c(76:126))

nir_14 <- read_stars(here("imagery", "subset_images", "rfl_2022_05_29.tif")) %>%
  st_set_dimensions(3, values = seq(1:425)) %>%
  filter(band %in% c(76:126)) %>%
  split(3) %>%
  setNames(c(76:126))

nir_all <- c(nir_1, nir_2, nir_3, nir_4, nir_5, nir_6, nir_7, nir_8, nir_9, nir_10, nir_11, nir_13, nir_14, along = 3) %>%
  st_set_dimensions(3, values = lubridate::ymd(c("2022_02_24", "2022_02_28", "2022_03_08", "2022_03_16", "2022_03_22", "2022_04_05", "2022_04_12", "2022_04_20", "2022_04_29", "2022_05_03", "2022_05_11", "2022_05_17", "2022_05_29")), names = "date") %>% 
  st_warp(dest = crsi_all) 

Adjusting elevation data

mllw_all <- mllw_all %>% 
  st_warp(dest = crsi_all) %>% 
  setNames("elevation")

Extracting Values based on week

soils_extract <- data.frame()

for (i in unique(lubridate::week(soils$paired_flight))){
  ndvi_week <- ndvi_all %>% 
    filter(lubridate::week(date) == i)
  
  mari_week <- mari_all %>% 
    filter(lubridate::week(date) == i)

  vssi_week <- vssi_all %>% 
    filter(lubridate::week(date) == i)
  
  crsi_week <- crsi_all %>% 
    filter(lubridate::week(date) == i)
  
  nir_week <- nir_all %>% 
    filter(lubridate::week(date) == i) %>% 
    st_as_sf
  
  soil_week <- soils %>% 
    mutate(week = lubridate::week(paired_flight)) %>% 
    filter(week == i)
  
  ndvi_extract <- ndvi_week %>% 
    st_extract(soil_week)
  
  mari_extract <- mari_week %>% 
    st_extract(soil_week)
  
  vssi_extract <- vssi_week %>% 
    st_extract(soil_week)
  
  crsi_extract <- crsi_week %>% 
    st_extract(soil_week)
  
  nir_extract <- soil_week %>% 
    st_join(nir_week) %>% 
    select("76":"126") %>% 
    st_drop_geometry()
  
  soils_week <- bind_cols(soil_week, ndvi_extract$ndvi, mari_extract$mari, vssi_extract$vssi, crsi_extract$crsi, nir_extract) %>% 
    rename("ndvi" = "...11",
           "mari" = "...12",
           "vssi" = "...13",
           "crsi" = "...14")
    
  soils_extract <- rbind(soils_extract, soils_week)
}

PCA

PCA Object

# Read in the data
nir_pcd <- soils_extract %>% # call object
  as.data.frame() %>% 
  select(electro_cond_mS_per_cm,
         elevation,
         ndvi:crsi,
         `76`:`126`) %>% # select the variables we want to explore
  drop_na()  # drops columns with NAs as PCA only works with numeric values

Perform the PCA

# Performing the PCA
nir_pca <- nir_pcd %>% # pipe in the data
  prcomp(scale. = TRUE) # rescaling of data is performed

Summarize

summary(nir_pca)

## Importance of components:
##                           PC1     PC2     PC3     PC4     PC5     PC6     PC7
## Standard deviation     7.1759 1.75884 1.06655 0.77939 0.64010 0.36867 0.22789
## Proportion of Variance 0.9034 0.05427 0.01996 0.01066 0.00719 0.00238 0.00091
## Cumulative Proportion  0.9034 0.95767 0.97763 0.98829 0.99548 0.99786 0.99877
##                            PC8     PC9    PC10    PC11    PC12    PC13    PC14
## Standard deviation     0.17019 0.08872 0.07819 0.06710 0.05182 0.04731 0.04508
## Proportion of Variance 0.00051 0.00014 0.00011 0.00008 0.00005 0.00004 0.00004
## Cumulative Proportion  0.99928 0.99942 0.99952 0.99960 0.99965 0.99969 0.99973
##                           PC15    PC16    PC17    PC18    PC19    PC20    PC21
## Standard deviation     0.03938 0.03687 0.03189 0.03122 0.02883 0.02796 0.02572
## Proportion of Variance 0.00003 0.00002 0.00002 0.00002 0.00001 0.00001 0.00001
## Cumulative Proportion  0.99975 0.99978 0.99979 0.99981 0.99983 0.99984 0.99985
##                           PC22    PC23    PC24    PC25    PC26    PC27    PC28
## Standard deviation     0.02501 0.02379 0.02299 0.02248 0.02156 0.02062 0.02012
## Proportion of Variance 0.00001 0.00001 0.00001 0.00001 0.00001 0.00001 0.00001
## Cumulative Proportion  0.99986 0.99987 0.99988 0.99989 0.99990 0.99991 0.99991
##                           PC29    PC30    PC31    PC32    PC33   PC34    PC35
## Standard deviation     0.01975 0.01867 0.01835 0.01813 0.01700 0.0165 0.01604
## Proportion of Variance 0.00001 0.00001 0.00001 0.00001 0.00001 0.0000 0.00000
## Cumulative Proportion  0.99992 0.99993 0.99993 0.99994 0.99994 1.0000 0.99995
##                           PC36    PC37    PC38    PC39    PC40    PC41   PC42
## Standard deviation     0.01574 0.01485 0.01443 0.01417 0.01311 0.01259 0.0125
## Proportion of Variance 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.0000
## Cumulative Proportion  0.99996 0.99996 0.99996 0.99997 0.99997 0.99997 1.0000
##                           PC43    PC44    PC45    PC46    PC47    PC48     PC49
## Standard deviation     0.01182 0.01164 0.01121 0.01095 0.01063 0.01046 0.009629
## Proportion of Variance 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.000000
## Cumulative Proportion  0.99998 0.99998 0.99998 0.99999 0.99999 0.99999 0.999990
##                            PC50     PC51     PC52    PC53     PC54     PC55
## Standard deviation     0.009267 0.009191 0.008831 0.00823 0.007837 0.006979
## Proportion of Variance 0.000000 0.000000 0.000000 0.00000 0.000000 0.000000
## Cumulative Proportion  0.999990 0.999990 1.000000 1.00000 1.000000 1.000000
##                            PC56     PC57
## Standard deviation     0.006765 0.006138
## Proportion of Variance 0.000000 0.000000
## Cumulative Proportion  1.000000 1.000000

Generate Biplot

#Creating the biplot via autoplot()

autoplot(nir_pca, # pca data
         x = 1,
         y = 2,
         data = nir_pcd, # original data
         loadings = TRUE, # show loadings
         loadings.label = TRUE, # label loadings
         loadings.colour = "black", # loading color
         loadings.label.colour = "black", # loading label color
         loadings.label.vjust = -.5) + # vertical justification
  ggtitle("Principle Component 1 and 2 Biplot") + # title
  labs(x = "Principle Component 1", 
       y= "Principle Component 2")+ # axis labels 
   theme(plot.title = element_text(color = "#5b4f41", size = 16),
            plot.background = element_rect("white"),
            panel.background = element_rect("#faf7f2"),
            panel.grid = element_line(linetype= "longdash",
                                      color = "#f0ece1"),
            axis.text = element_text(color = "#5b4f41"),
            axis.title = element_text(color = "#5b4f41"),
            strip.background = element_rect("white"),
            axis.line = element_line(color = "#5b4f41")) # change themes

autoplot(nir_pca, # pca data
         x = 1,
         y = 3,
         data = nir_pcd, # original data
         loadings = TRUE, # show loadings
         loadings.label = TRUE, # label loadings
         loadings.colour = "black", # loading color
         loadings.label.colour = "black", # loading label color
         loadings.label.vjust = -.5) + # vertical justification
  ggtitle("Principle Component 1 and 3 Biplot") + # title
  labs(x = "Principle Component 1", 
       y= "Principle Component 3")+ # axis labels 
   theme(plot.title = element_text(color = "#5b4f41", size = 16),
            plot.background = element_rect("white"),
            panel.background = element_rect("#faf7f2"),
            panel.grid = element_line(linetype= "longdash",
                                      color = "#f0ece1"),
            axis.text = element_text(color = "#5b4f41"),
            axis.title = element_text(color = "#5b4f41"),
            strip.background = element_rect("white"),
            axis.line = element_line(color = "#5b4f41")) # change themes

autoplot(nir_pca, # pca data
         x = 1,
         y = 4,
         data = nir_pcd, # original data
         loadings = TRUE, # show loadings
         loadings.label = TRUE, # label loadings
         loadings.colour = "black", # loading color
         loadings.label.colour = "black", # loading label color
         loadings.label.vjust = -.5) + # vertical justification
  ggtitle("Principle Component 1 and 4 Biplot") + # title
  labs(x = "Principle Component 1", 
       y= "Principle Component 4")+ # axis labels 
   theme(plot.title = element_text(color = "#5b4f41", size = 16),
            plot.background = element_rect("white"),
            panel.background = element_rect("#faf7f2"),
            panel.grid = element_line(linetype= "longdash",
                                      color = "#f0ece1"),
            axis.text = element_text(color = "#5b4f41"),
            axis.title = element_text(color = "#5b4f41"),
            strip.background = element_rect("white"),
            axis.line = element_line(color = "#5b4f41")) # change themes

autoplot(nir_pca, # pca data
         x = 2,
         y = 3,
         data = nir_pcd, # original data
         loadings = TRUE, # show loadings
         loadings.label = TRUE, # label loadings
         loadings.colour = "black", # loading color
         loadings.label.colour = "black", # loading label color
         loadings.label.vjust = -.5) + # vertical justification
  ggtitle("Principle Component 2 and 3 Biplot") + # title
  labs(x = "Principle Component 2", 
       y= "Principle Component 3")+ # axis labels 
   theme(plot.title = element_text(color = "#5b4f41", size = 16),
            plot.background = element_rect("white"),
            panel.background = element_rect("#faf7f2"),
            panel.grid = element_line(linetype= "longdash",
                                      color = "#f0ece1"),
            axis.text = element_text(color = "#5b4f41"),
            axis.title = element_text(color = "#5b4f41"),
            strip.background = element_rect("white"),
            axis.line = element_line(color = "#5b4f41")) # change themes

autoplot(nir_pca, # pca data
         x = 2,
         y = 4,
         data = nir_pcd, # original data
         loadings = TRUE, # show loadings
         loadings.label = TRUE, # label loadings
         loadings.colour = "black", # loading color
         loadings.label.colour = "black", # loading label color
         loadings.label.vjust = -.5) + # vertical justification
  ggtitle("Principle Component 2 and 4 Biplot") + # title
  labs(x = "Principle Component 2", 
       y= "Principle Component 4")+ # axis labels 
   theme(plot.title = element_text(color = "#5b4f41", size = 16),
            plot.background = element_rect("white"),
            panel.background = element_rect("#faf7f2"),
            panel.grid = element_line(linetype= "longdash",
                                      color = "#f0ece1"),
            axis.text = element_text(color = "#5b4f41"),
            axis.title = element_text(color = "#5b4f41"),
            strip.background = element_rect("white"),
            axis.line = element_line(color = "#5b4f41")) # change themes

autoplot(nir_pca, # pca data
         x = 3,
         y = 4,
         data = nir_pcd, # original data
         loadings = TRUE, # show loadings
         loadings.label = TRUE, # label loadings
         loadings.colour = "black", # loading color
         loadings.label.colour = "black", # loading label color
         loadings.label.vjust = -.5) + # vertical justification
  ggtitle("Principle Component 3 and 4 Biplot") + # title
  labs(x = "Principle Component 3", 
       y= "Principle Component 4")+ # axis labels 
   theme(plot.title = element_text(color = "#5b4f41", size = 16),
            plot.background = element_rect("white"),
            panel.background = element_rect("#faf7f2"),
            panel.grid = element_line(linetype= "longdash",
                                      color = "#f0ece1"),
            axis.text = element_text(color = "#5b4f41"),
            axis.title = element_text(color = "#5b4f41"),
            strip.background = element_rect("white"),
            axis.line = element_line(color = "#5b4f41")) # change themes

# Creating screeplot
sd_vec <- nir_pca$sdev # standard deviation
var_vec <- sd_vec^2 # variance

pc_names <- colnames(nir_pca$rotation) # names of PCs

pct_expl_df <- data.frame(v = var_vec, # new data frame for screeplot; variance
                          pct_v = var_vec / sum(var_vec), # percent of variance
                          pc = fct_inorder(pc_names)) %>%  # orders rows
  mutate(pct_label = paste0(round(pct_v * 100, 1), "%")) # adds character %

ggplot(pct_expl_df, aes(x = pc, y = v))+ # graphs via PC and variance
  geom_col()+
  geom_text(aes(label = pct_label), vjust = 0, nudge_y = 0.008)+
  labs(x = "Principle Components",
       y = "Variances") +
  ggtitle("PCA Screeplot")+
  theme(plot.title = element_text(hjust = .5))

Fitting model with all variables currently available

Fitting the model to all variables

fit_all <- train(electro_cond_mS_per_cm ~ elevation + ndvi + mari + vssi + crsi, 
                 data = soils_extract, 
                 method = "rf", 
                 trControl = trainControl(method = "cv", number = 10), 
                 importance = TRUE)

Evaluating Model

# using caret
fit_all_imp <- as.data.frame(fit_all$finalModel$importance)
fit_all_imp_scaled <- predict(preProcess(fit_all_imp, method = c("range")), fit_all_imp) %>% 
  rownames_to_column(var = "variable")

ggplot(fit_all_imp_scaled) +
  geom_point(aes(x = IncNodePurity, y = variable), color = "#E69512")+
  geom_segment(aes(y= variable, x = 0, yend = variable, xend = IncNodePurity), color = "#E69512", size = 1.2)+
  geom_point(aes(x = fit_all_imp_scaled$`%IncMSE`, y = variable), color ="#3A5D3D")+ 
  geom_segment(aes(y= variable, x = 0, yend = variable, xend = fit_all_imp_scaled$`%IncMSE`), color = "#3A5D3D", linetype = "dashed")+
   labs(y = "Variables",#labels
       x = "Importance (Scaled)")+
  ggtitle("Random Forest Variable Importance")+#plot title
  theme(plot.title = element_text(color = "#5b4f41", size = 16, hjust = 0.5),#plot themes
        plot.background = element_rect("white"),
        panel.background = element_rect("#faf7f2"),
        panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
        axis.text = element_text(color = "#5b4f41"),
        axis.title = element_text(color = "#5b4f41"),
        strip.background = element_rect("#f8f8f8"),
        strip.text = element_text(color = "#5b4f41"),
        axis.line = element_line(color = "#5b4f41"))

cor(fit_all$finalModel$y, fit_all$finalModel$predicted)

## [1] 0.8421477

Ploting Model Fit

ggplot()+#ggplot object
  geom_abline(slope = 1, intercept = 0, linetype = "dotted")+
  geom_point(aes(y = fit_all$finalModel$y, x = fit_all$finalModel$predicted), color = "#E69512")+#data to make points from
  geom_smooth(method = "lm", aes(y= fit_all$finalModel$y, x = fit_all$finalModel$predicted), color =  "#3A5D3D")+#smooth line data
  scale_x_continuous(expand = c(0, 0), limits = c(0, 14)) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 14))+
  labs(y = "Actual Electrical Conductivity (mS/cm)",#labels
       x = "Predicted Electrical conductivity (mS/cm)")+
  ggtitle("Elevation and Indices (caret method)")+#plot title
  theme(plot.title = element_text(color = "#5b4f41", size = 16, hjust = 0.5),#plot themes
        plot.background = element_rect("white"),
        panel.background = element_rect("#faf7f2"),
        panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
        axis.text = element_text(color = "#5b4f41"),
        axis.title = element_text(color = "#5b4f41"),
        strip.background = element_rect("#f8f8f8"),
        strip.text = element_text(color = "#5b4f41"),
        axis.line = element_line(color = "#5b4f41"))

#ggsave(filename = "all_plot_caret.png", plot = last_plot(), width = 12, height = 8, device = "png", dpi = 500) #write an image of the plot

all_corr <- cor(y= fit_all$finalModel$y, x = fit_all$finalModel$predicted)
all_corr

## [1] 0.8421477

Fitting the model to spectral Indices

fit_si <- train(electro_cond_mS_per_cm ~ ndvi + mari + vssi + crsi + `76` + `77` + `78` + `79` + `80` + `81` + `82` + `83` + `84` + `85` + `86` + `87` + `88` + `89` + `90` + `91` + `92` + `93` + `94` + `95` + `96` + `97` + `98` + `99` + `100` + `101` + `102` + `103` + `104` + `105` + `106` + `107` + `108` + `109` + `110` + `111` + `112` + `113` + `114` + `115` + `116` + `117` + `118` + `119` + `120` + `121` + `122` + `123` + `124` + `125` + `126`, 
                 data = soils_extract, 
                 method = "rf", 
                 trControl = trainControl(method = "cv", number = 10), 
                 importance = TRUE)

Evaluating Model

fit_si_imp <- as.data.frame(fit_si$finalModel$importance)
fit_si_imp_scaled <- predict(preProcess(fit_si_imp, method = c("range")), fit_si_imp) %>% 
  rownames_to_column(var = "variable")


ggplot(fit_si_imp_scaled) +
  geom_point(aes(x = IncNodePurity, y = variable), color = "#E69512")+
  geom_segment(aes(y= variable, x = 0, yend = variable, xend = IncNodePurity), color = "#E69512", size = 1.2)+
  geom_point(aes(x = fit_si_imp_scaled$`%IncMSE`, y = variable), color ="#3A5D3D")+ 
  geom_segment(aes(y= variable, x = 0, yend = variable, xend = fit_si_imp_scaled$`%IncMSE`), color = "#3A5D3D", linetype = "dashed")+
   labs(y = "Variables",#labels
       x = "Importance (Scaled)")+
  ggtitle("Random Forest Variable Importance")+#plot title
  theme(plot.title = element_text(color = "#5b4f41", size = 16, hjust = 0.5),#plot themes
        plot.background = element_rect("white"),
        panel.background = element_rect("#faf7f2"),
        panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
        axis.text = element_text(color = "#5b4f41"),
        axis.title = element_text(color = "#5b4f41"),
        strip.background = element_rect("#f8f8f8"),
        strip.text = element_text(color = "#5b4f41"),
        axis.line = element_line(color = "#5b4f41"))

cor(y = fit_si$finalModel$y, x= fit_si$finalModel$predicted)

## [1] 0.5060143

Observed v Predicted Plot

ggplot()+#ggplot object
  geom_abline(slope = 1, intercept = 0, linetype = "dotted")+
  geom_point(aes(y= fit_si$finalModel$y, x = fit_si$finalModel$predicted), color = "#E69512")+#data to make points from
  geom_smooth(method = "lm", aes(y= fit_si$finalModel$y, x = fit_si$finalModel$predicted), color =  "#3A5D3D")+#smooth line data
  scale_x_continuous(expand = c(0, 0), limits = c(0, 14)) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 14))+
  labs(y = "Actual Electrical Conductivity (mS/cm)",#labels
       x = "Predicted Electrical conductivity (mS/cm)")+
  ggtitle("Spectral Indices")+#plot title
  theme(plot.title = element_text(color = "#5b4f41", size = 16, hjust = 0.5),#plot themes
        plot.background = element_rect("white"),
        panel.background = element_rect("#faf7f2"),
        panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
        axis.text = element_text(color = "#5b4f41"),
        axis.title = element_text(color = "#5b4f41"),
        strip.background = element_rect("#f8f8f8"),
        strip.text = element_text(color = "#5b4f41"),
        axis.line = element_line(color = "#5b4f41"))

#ggsave(filename = "si_plot.png", plot = last_plot(), width = 12, height = 8, device = "png", dpi = 500) #write an image of the plot

Random Forest regression of Soil salinity and elevation

fit_ele <- train(electro_cond_mS_per_cm ~ elevation, data = soils, method = "rf", trControl = trainControl(method = "cv", number = 10), importance = T)#random forest fitting of the data with stated equation

Correlation

elevation_cor <- cor(y = fit_ele$finalModel$y, x = fit_ele$finalModel$predicted)#calculate correlation coefficient
elevation_cor#call/print correlation value

## [1] 0.9208296

Observed v Predicted

ggplot()+#ggplot object
  geom_abline(slope = 1, intercept = 0, linetype = "dotted")+
  geom_point(aes(y = fit_ele$finalModel$y, x = fit_ele$finalModel$predicted), color = "#E69512")+#data to make points from
  geom_smooth(method = "lm", aes(y = fit_ele$finalModel$y, x = fit_ele$finalModel$predicted), color =  "#3A5D3D")+#smooth line data
  scale_x_continuous(expand = c(0, 0), limits = c(0, 14)) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 14))+
  labs(y = "Actual Electrical Conductivity (mS/cm)",#labels
       x = "Predicted Electrical conductivity (mS/cm)")+
  ggtitle("Elevation")+#plot title
  theme(plot.title = element_text(color = "#5b4f41", size = 16, hjust = 0.5),#plot themes
        plot.background = element_rect("white"),
        panel.background = element_rect("#faf7f2"),
        panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
        axis.text = element_text(color = "#5b4f41"),
        axis.title = element_text(color = "#5b4f41"),
        strip.background = element_rect("#f8f8f8"),
        strip.text = element_text(color = "#5b4f41"),
        axis.line = element_line(color = "#5b4f41"))

#ggsave(filename = "ele_plot.png", plot = last_plot(), width = 12, height = 8, device = "png", dpi = 500) #write an image of the plot

Compiling all raster data

raster_all <- c(crsi_all, ndvi_all, vssi_all, mari_all, mllw_all)

Predicting values

all_pred <- predict(raster_all, fit_all, drop_dimensions = F)

Making a Map

all_pred_feb <- all_pred %>% 
  filter(date == "2022-02-24") %>% 
  st_crop(y = boundary, crop = TRUE, as_points = FALSE)


ggplot()+
  geom_stars(data = all_pred_feb, aes(x, y, fill = prediction)) +
  scale_fill_gradientn(colors = rev(calecopal::cal_palette("desert", 2, "continuous")))+
  labs(y = "Northing",#labels
       x = "Easting",
       fill = "Electrical 
Conductivity 
(mS/cm)")+
  ggtitle("Soil Salinity")+#plot title
  theme(plot.title = element_text(color = "#5b4f41", size = 16, hjust = 0.5),#plot themes
        plot.background = element_rect("white"),
        panel.background = element_rect("#faf7f2"),
        panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
        axis.text = element_text(color = "#5b4f41"),
        axis.title = element_text(color = "#5b4f41"),
        strip.background = element_rect("#f8f8f8"),
        strip.text = element_text(color = "#5b4f41"),
        axis.line = element_line(color = "#5b4f41"),
        legend.title = element_text(color = "#5b4f41"),
        legend.text = element_text(color = "#5b4f41"))

# Interactive Map

mapview::mapview(all_pred_feb)

TESTING OF PLSR INSTEAD!

Fitting the model to all variables

fit_all_pls <- train(electro_cond_mS_per_cm ~ elevation + ndvi + mari + vssi + crsi, 
                 data = soils_extract, 
                 method = "pls", 
                 trControl = trainControl(method = "cv", number = 10), 
                 importance = TRUE)

Evaluating Model

# using caret

predicted_values <- as.data.frame(fit_all_pls$finalModel$fitted.values) %>% 
  select(`.outcome.3 comps`)

training_data <- as.data.frame(fit_all_pls$trainingData$.outcome)

train_pred <- c(predicted_values, training_data) %>% 
  as.data.frame()

cor(training_data, predicted_values)

##                                   .outcome.3 comps
## fit_all_pls$trainingData$.outcome        0.6955281

Ploting Model Fit

ggplot()+#ggplot object
  geom_abline(slope = 1, intercept = 0, linetype = "dotted")+
  geom_point(data = train_pred, aes(y = fit_all_pls.trainingData..outcome, x = `.outcome.3.comps`), color = "#E69512")+#data to make points from
  geom_smooth(method = "lm", data = train_pred, aes(y = fit_all_pls.trainingData..outcome, x = `.outcome.3.comps`), color =  "#3A5D3D")+#smooth line data
  labs(y = "Actual Electrical Conductivity (mS/cm)",#labels
       x = "Predicted Electrical conductivity (mS/cm)")+
  ggtitle("Elevation and Indices (PLSR)")+#plot title
  theme(plot.title = element_text(color = "#5b4f41", size = 16, hjust = 0.5),#plot themes
        plot.background = element_rect("white"),
        panel.background = element_rect("#faf7f2"),
        panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
        axis.text = element_text(color = "#5b4f41"),
        axis.title = element_text(color = "#5b4f41"),
        strip.background = element_rect("#f8f8f8"),
        strip.text = element_text(color = "#5b4f41"),
        axis.line = element_line(color = "#5b4f41"))

#ggsave(filename = "all_plot_caret.png", plot = last_plot(), width = 12, height = 8, device = "png", dpi = 500) #write an image of the plot

all_pred_pls <- predict(raster_all, fit_all_pls, drop_dimensions = F)

all_pred_feb_pls <- all_pred_pls %>% 
  filter(date == "2022-02-24") %>% 
  st_crop(y = boundary, crop = TRUE, as_points = FALSE)


ggplot()+
  geom_stars(data = all_pred_feb_pls, aes(x, y, fill = prediction)) +
  scale_fill_gradientn(colors = rev(calecopal::cal_palette("desert", 2, "continuous")))+
  labs(y = "Northing",#labels
       x = "Easting",
       fill = "Electrical 
Conductivity 
(mS/cm)")+
  ggtitle("Soil Salinity")+#plot title
  theme(plot.title = element_text(color = "#5b4f41", size = 16, hjust = 0.5),#plot themes
        plot.background = element_rect("white"),
        panel.background = element_rect("#faf7f2"),
        panel.grid = element_line(linetype= "longdash", color = "#f0ece1"),
        axis.text = element_text(color = "#5b4f41"),
        axis.title = element_text(color = "#5b4f41"),
        strip.background = element_rect("#f8f8f8"),
        strip.text = element_text(color = "#5b4f41"),
        axis.line = element_line(color = "#5b4f41"),
        legend.title = element_text(color = "#5b4f41"),
        legend.text = element_text(color = "#5b4f41"))

Interactive Map

mapview::mapview(all_pred_feb_pls)